Cosmological N-Body simulation: Techniques, Scope and Status 
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Cosmological N-Body simulations have become an essential tool for studying formation of large 
scale structure. These simulations are computationally challenging even though the available com- 
puting power gets better every year. A number of efficient algorithms have been developed to run 
large simulations with better dynamic range and resolution. We discuss key algorithms in this 
review, focusing on techniques used and their efficacy. N-Body simulations solve a model that is 
an approximation of the physical model to be simulated, we discuss limitations arising from this 
approximation and techniques employed for solving equations. Apart from simulating models of 
structure formation, N-Body simulations have also been used to study aspects of gravitational clus- 
tering. Simulating formation of galaxies requires us to take many physical process into account; we 
review numerical implementations of key processes. 



I. INTRODUCTION 

Large scale structures like galaxies and clusters of 
galaxies are believed to have formed by amplification of 
small perturbations |Jj,i^iSfc«4n^«^0iBl3- Galaxies are 
highly over-dense systems, matter density p in galaxies 
is thousands of times larger than the average density p 
in the universe. Typical density contrast {5 = p/p—l) 
in matter at these scales in the early universe, e.g. at 
the time of decoupling of matter and radiation was much 
smaller than unity. Thus the problem of galaxy formation 
and the large scale distribution of galaxies is essentially 
one of evolving density perturbations from small initial 
values to the large values we encounter today. 

The universe is assumed to be homogeneous and 
isotropic at large scales and this assumption is consis- 
tent with observations and there are strong limits on 
departures from homogeneity and isotropy A ho- 

mogeneous and isotropic universe is described by Fried- 
man equations in the general theory of relativity. As 
long as perturbations in the gravitational potential are 
small, we can treat density fluctuations as perturbations 
about a Friedman universe. If perturbations are in non- 
relativistic matter, as appears to be the case in our uni- 
verse, we can work in the Newtonian limit for studying 
their evolution. The back reaction of perturbations on 
the universe is not taken into account, i.e., we do not 
worry about the effect of perturbations on the average 
global properties of the universe. Studies have shown 
that this back- reaction is ignorable in most cases [TllIT^ . 
Local variations in density, etc. can lead to dispersion in 
values of cosmological parameters determined through lo- 
cal observations 13] but this problem can be controlled 
by making measurements at larger scales. 

Gravity is the dominant force at large scales and is be- 
lieved to drive growth of perturbations. Magnetic field 
is the only other interaction that can lead to formation 
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of large scale structures, this leads to very distinct sig- 
natures and does not appear to be the dominant factor 
in our universe 0, 0, ■ We will assume that gravity 
is the only interaction responsible for growth of pertur- 
bations at large scales. Equations that describe the evo- 
lution of density perturbations in non-relativistic matter 
due to gravitational interaction in an expanding universe 
have been known for a long time [T^ . These equations 
can be solved analytically for small density contrasts, and 
for highly symmetric situations. But apart from such spe- 
cial cases, few solutions are known. Many approximate 
solutions arc known ^IS„ia,^^, 22, 23, 2424| and 
are useful in understanding the evolution of perturba- 
tions in the quasi-linear regime, these fail when density 
contrast becomes large [5 3> 1). Cosmological N-Body 
simulations are an essential tool for evolving density per- 
turbations in the non-linear regime. Fluctuations in the 
gravitational potential do not grow by a large amount 
even as density contrast increases by several orders of 
magnitude [2^l23| . therefore the Newtonian approxima- 
tion continues to be a valid framework. At galactic scales, 
gas dynamical and other effects play an important role 
and need to be taken into account for a detailed solution 
of the problem. 

N-Body simulations are used for a variety of applica- 
tions. Simulations of specific models of dark matter allow 
us to make predictions for these models and compare with 
observations. Simulations allow us to carry out numerical 
experiments with initial conditions that have little to do 
with the real universe. The purpose of such experiments 
is to understand the physics of gravitational collapse in 
an expanding universe. Simulations are also used for test- 
ing approximate solutions for growth of density pertur- 
bations, comparisons with N-Body simulations allow us 
to validate these approximations and understand when 
these approximations are useful. Lastly, we can calibrate 
methods for analysing observations on mock catalogues 
made from N-Body simulations. We can test whether a 
particular method works or not because in N-Body all 
the details are known whereas the same is not true of the 
real universe. 

The physical parameters of the problem make cosmo- 



2 



logical simulations a challenging task. Unlike simulations 
of systems like globular clusters that can be treated as 
(relatively) isolated systems, the universe does not have a 
boundary. Therefore cosmological simulations need to be 
run with periodic boundary conditions. An exception are 
simulations of large spherical volumes that do not suffer 
significant deformation during the course of evolution. 

Observations suggest that density perturbations are 
present at all scales that have been probed by observa- 
tions. Amplitude of fluctuations at small scales is large 
and it drops at large scales (/ ^ 10 Mpc where 1 Mpc = 
3.08 X 10^** cm.). Figure n] shows the root mean square 
fluctuations in mass as a function of length scale for one of 
the popular models. Fluctuations at scales much larger 
than 100 Mpc are generally not relevant for structure 
formation at the scale of galaxies. Thus the physical size 
corresponding to the periodic simulation box should be 
at least 100 Mpc, unless these are meant for studying 
large scale structure at early times when the amplitude 
of fluctuations is small and a smaller simulation box is 
acceptable. We will revisit this issue and discuss it more 
quantitatively but this approximate figure will sufhcc at 
present. 

If we wish to study the distribution of galaxies in detail 
then the mass of each particle in the simulation should 
be much smaller than the mass of a typical galaxy. This, 
with the requirement that the simulation box should be 
more than 100 Mpc across implies that the N-Body sim- 
ulation must be done with at least 10^ particles. The 
large number of particles required is one of the things 
that makes cosmological simulations challenging. 

Observations show that the dominant component of 
matter in the universe does not radiate light and thus 
cannot be seen, except through its gravitational effect 
on visible matter Observational evidence points 

towards non-relativistic dark matter [27j .2^ J9J. This 
makes our task simple as motions of non-relativistic mat- 
ter can be studied in the fluid limit at large scales. 
(Relativistic dark matter will have significant pres- 
sure/velocity dispersion and one cannot take the fluid 
limit. In such a case one needs to solve the Boltzmann 
equation in order to correctly model the effects of free 
streaming Is^l-) Little is known about what constitutes 
this dark matter though it is often assumed to be made of 
Weakly Interacting Massive Particles (WIMP). There are 
strong limits on the interaction cross section of WIMPs 
from astrophysical considerations [s^l- Non-relativistic, 
non-interacting (coUisionless) dark matter is known as 
Cold Dark Matter (CDM) [3l|3l|3|. Dark matter in- 
teracts only through gravity and that makes simulating 
growth of perturbations somewhat simpler. 

Dark matter dominates over normal matter in terms 
of density by a large factor 28] , therefore N-Body sim- 
ulations with the entire matter density in dark matter 
provide a good first approximation for the distribution 
of matter. In any case, gravity is the dominant interac- 
tion and normal matter is expected to follow dark matter 
at large scales. 



Cosmological simulations differ from other types of N- 
Body simulations as the background is not static and ex- 
pansion of the universe has to be taken into account. It is 
convenient to work in comoving coordinates that expand 
with the universe. Equations in comoving coordinates 
deal with the evolution of perturbations and the average 
quantities (density, velocity, etc.) are scaled out. We will 
discuss this in greater detail in the following section. 

In N-Body simulations, each particle represents a very 
large number of dark matter particles and interaction of 
two particles in N-Body simulations should mimic the 
interaction of two "fluid elements" . The fluid elements 
being simulated have a physical size and at scales com- 
parable to this, the fluid elements should feel much less 
force than two point particles. This is done by assuming 
that the particles have a finite size and density profile, 
this leads to an effective softening of force at small scales. 
Clearly, the force should be softened at scales compara- 
ble to the (local) average inter-particle separation in the 
N-Body simulation. A much smaller softe ning length can 
lead to unwanted two body relaxation |35ll36| |. The form 
of force softening also plays an important role, force soft- 
ened in such a way that it matches inverse square force 
beyond the softening length is better [s^,!!^- If the soft- 
ened force approaches inverse square force only asymp- 
totically then the difference is equivalent to error in force 
at large separations and can lead to spurious effects. 

This ends the overview of the physical requirements 
and associated approximations for cosmological simula- 
tions. The next section contains a detailed discussion 
of cosmological N-Body simulations that take only grav- 
itational interaction into account. It is very important 
to understand limitations of N-Body simulations as us- 
ing simulations outside the domain of validity can easily 
lead to incorrect results. Limitations of N-Body simu- 
lations are discussed in the following section. N-Body 
simulations have been used to further our understand- 
ing of gravitational clustering in the non-linear regime, 
we discuss some relevant issues. This is followed by an 
overview of other processes that must be taken into ac- 
count for a complete study of galaxy formation, we also 
discuss numerical implementations of these physical pro- 
cesses in N-Body simulations. 



II. GRAVITY ONLY SIMULATIONS 

In this section we discuss N-Body simulations where 
only gravitational interaction is taken into account. In 
an N-Body simulation a given density- velocity field is rep- 
resented by a set of particles [S^- Density as a function 
of position is obtained by averaging over this distribution 
of particles. 

Evaluation of force and solving the equation of motion 
are two key components of N-Body simulations. Setting 
up relevant initial conditions for cosmological simulations 
is another important aspect. We discuss the equation of 
motion and its integration in the following subsection. 
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This is followed by a discussion of algorithms for force 
calculations. We end this section with a discussion on 
setting up the initial conditions. 

A. Equation of Motion 

The equations that govern the evolution of a given 
distribution of particles are obtained by starting with 
the standard (Newtonian) equation of motion for a 
set of gravitationally interacting particles in the physi- 
cal/proper coordinates. Transforming to comoving co- 
ordinates allows us to focus on perturbations in density 
and velocity. Information about expansion of the uni- 
verse appears in these equations through the scale factor 
aUy, obtained as a solution to the Friedman equations 
[llalllllSllllll- The scale factor and other cosmo- 
logical parameters in these equations carry information 
about the average quantities that are scaled out in the 
coordinate transformation. 

x + = -4tV(/) (1) 
a 

= 'iTTGp(t)a'^S = -H^no- (2) 
2 a 

Here x is the comoving position of a particle and is re- 
lated to the physical/proper position r = a(t)x, with a{t) 
being the scale factor, (j) is the gravitational potential 
due to density perturbations, Hq is the Hubble's con- 
stant and Qq is the density parameter of non-relativistic 
matter at the present epoch ^L, 2, 3.^ M- 
sumed that the relativistic components do not cluster, or 
are negligible in our universe. These equations are valid 
for non-relativistic matter (w <C c, (/> <C c^) at scales that 
are much smaller than the Hubble radius (/ ^ c/Hq) 

m 

The expansion of the universe acts as a viscous force 
in comoving coordinates. This drag opposes gravitational 
infall and as a result the growth of density perturbations 
is slower in an expanding universe. The time scale over 
which gravitational infall occurs (in absence of expan- 
sion) is comparable with the expansion time scale (a/a), 
therefore velocities of particles do not become very large 
during infall. As a result integrating equation of motion 
is simpler in cosmological N-Body simulations. 

The basic idea for numerical integration is as follows. 
The equation of motion expresses the second derivative 
of position in terms of position, velocity and time. Posi- 
tion and velocity at later times are expressed in terms of 
position and velocity at earlier times using a truncated 
Taylor series. The simplest truncation is not sufficiently 
accurate and the resulting error is of order in one 
time step, where h is the time step |3^0,0. The key 
constraint in cosmological simulations is that force evalu- 
ation is very time consuming and one wishes to minimise 
the number of force evaluations per time step. Mainly 
due to this reason, cosm olog ical N-Body simulations use 
the Leap-Frog method |39l liol l4ll | for integrating the 



equation of motion as it requires only one evaluation of 
force and the error is of order h^. Performance of the 
Leap-Frog integrator can be improved considerably by 
making small modifications 42], but such modifications 
are often more useful in non-cosmological N-Body simu- 
lations. 

Time step h is typically chosen so that momentum is 
conserved and energy evolves according to the Irvine- 
Layzer equation ^3)I3E1|- Monitoring consistency with 
the Irvine-Layzer equation requires care and adds signif- 
icantly to the number of operations to be carried out in 
an N-Body simulation , hence it is usual to carry out 
test runs and fix the value of h. Additional tests can be 
devised, e.g. we can require that the clustering of power 
law models evolves in a scale invariant manner even in 
the strongly non-linear regime. 

Optimum value of time step h depends on the dis- 
tribution of particles and it changes as this distribution 
evolves. It is common to use a time step that varies with 
time so that the N-Body code does not use too small a 
time step when a larger value will do, or use too large an 
h when a smaller value is required for conserving integrals 
of motion. It is possible to generalise even further and 
choose a different h for each particle as well, motivation 
for this being that a few particles in a very dense regions 
require a small h whereas most particles are not in such 
regions. There are several methods of implementing this 
in N-Body simulations, e.g. see Main considera- 

tion is to ensure that the positions and velocities of all 
particles are synchronised at frequent intervals. Using 
individual time steps can speed up N-Body simulations 
by a significant amount. 

B. Calculating Force 

Gravitational force in the Newtonian limit falls as 1 /r^, 
hence it is a long range force and we cannot ignore force 
due to distant particles. This makes calculation of force 
the most time consuming task in N-Body simulations. 
As a result, a lot of attention has been focused on this 
aspect and many algorithms and optimising schemes have 
been developed. We will discuss the major algorithms in 
some detail and briefiy summarise other developments. 
We refer the reader to ^3 for a detailed review. In the 
following discussion, we also review some algorithms that 
are not used in cosmological N-Body simulations as these 
can be used in hybrid algorithms. 

1. Direct Summation or Particle-Particle Method 

The most obvious approach to the problem of force 
calculation is to carry out a direct pairwise summation 
over all particles. This is also called the Particle-Particle 
(PP) method and this works very well for a small number 
of particles. Most early simulations used this method, 
see |43| for an early approach, 'sSj for a discussion of 
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characteristics of this method, and, [5l| for an overview 
of the direct summation method. Earhest cosmological 
simulations also used this method The number of 

terms in the pairwise summation increases in proportion 
with A^^, where N is the number of particles. This rapid 
variation limited the early cosmological simulations to 
about 10'^ particles. Limitations of the PP approach were 
realised as focus shifted to other methods. 

The last decade has seen a revival of sorts for this 
method with the advent of the GRAPE chip With 
the latest version of GRAPE ^] it is possible to simulate 
systems with more than 10^ particles. Development of ef- 
ficient parallel algorithms [53 has also made this method 
competitive. 

It is difficult to implement periodic boundary condi- 
tions in the PP method. The only optimised method 
available for this is the Ewald summation [s^. One of 
the first concrete proposals to use this method in cosmo- 
logical N-Body simulations is [s^l, though Ewald sum- 
mation had been used for testing other methods [4fil |. 
Adding periodic boundary conditions to a PP code re- 
mains a difficult proposition and this method is not used 
very often for cosmological N-Body simulations. 



2. Tree Method 

The key limitation of the PP method is the rapid in- 
crease in computational load with the number of parti- 
cles in the N-Body simulation. This in turn arises from 
adding individual contribution to the force due to each 
particle. The force of a distant group of particles can be 
approximated by the force due to a single pseudo parti- 
cle located at the centre of mass of the group, with mass 
equal to the total mass of the group of particles. This 
approximation changes the scaling of the number of cal- 
culations from N'^ to A^logiV. 

Efficient division of particles into groups can be done 
by arranging particles in a tree structure |58, 59J. The 
simulation volume is taken to be a cube and is divided 
into smaller cubes with 1/8 the volume each at every 
stage till the smallest cells have only one particle in them. 
Larger cells serve as groups of particles for a rapid calcu- 
lation of force. An essential ingredient is the criterion for 
deciding whether a group of particles can be considered 
distant or not. This is called the cell acceptance crite- 
rion and the error in approximation is controlled by the 
choice of this criterion lUii 13 See [HiO, El El 
for a detailed study of characteristics of the tree code, 
in particular of errors and timing as a function of the 
distribution of particles and the cell acceptance criteria. 

Accuracy of the tree approximation can be improved 
by retaining information about moments of the particle 
distribution in the group, e.g. the quadruple moment. 

The tree code can be optimised by vectorising the cal- 
culation of force [63l l63 | . The set of acceptable distant 
groups of particles is almost the same for neighbouring 
particles and sharing of interaction lists amongst neigh- 



bouring particles can further improve the performance 
of the tree code . The tree code can be parallelised 
efficiently ji^, and a parallel tree code has also been 
implemented using the GRAPE chip [63|. The parallel 
algorithm divides the simulation box into domains with 
equal number of particles and calculations for each do- 
main are done by a different processor. This scheme can 
be improved by dividing into domains with equal com- 
putational load j66i |. 

Periodic boundary conditions are difficult to imple- 
ment with tree codes, the level of difficulty being similar 
to that with the PP codes. Some innovative schemes have 
been tried |6^ . Implementations using the Ewald sum- 
mation [s^ add a large computational overhead ^3, . 
But in spite of these difficulties, tree codes have been used 
very effectively for cosmological N-Body simulations. 

3. Fast Multipole Method 

The performance of tree codes can be improved upon 
by using including higher moments of mass distribu- 
tion in cells. These and use of some other optimisation 
schemes leads to the fast multipole method. The number 
of computational operations in the fast multipole method 
[70| scales as A'', the number of particles. 

Inclusion of higher moments can be also modelled in 
terms of pseu do-p articles for easy implementation on the 
GRAPE chip llllll. 

An explicitly momentum conserving extension of the 
tree code has also been proposed and implemented 
^3', '7^ , here the number of computations required scales 
linearly with the number of particles. 

These codes also suffer from the problem of open 
boundary conditions and cannot be adapted very easily 
to cosmological problems. 

4. Particle-Mesh Method 

Particle-Mesh (PM) method jH El III S III S IH 

has been used extensively for cosmological simulations 
and was the first method to be used for "large" [N ~ 10^) 
simulations [z^. In PM codes, the fact that the Poisson's 
equation (eqn|51) is a simple algebraic equation in Fourier 
space is combined with Fast Fourier Transforms (EFT) 
I40II41I . FFT requires sampling of functions at uniformly 
spaced points, and a grid/mesh is used for this. Usually 
the simulation volume is taken to be a cube with equal 
number of grid/mesh points along each axis. 

Particles are used for representing the density and ve- 
locity field and we compute density at grid points by 
using weight functions 39]. Density contrast 5 and the 
potential <j) is defined on the grid for solving the Poisson 
equation. Use of particles and a mesh gives this method 
the very appropriate name. 

By using Fourier methods we get periodic boundary 
conditions for free. The use of a mesh also softens the 
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force naturally at small scales, though this softening leads 
to underestimation of force to fairly large scales |73, UM ■ 
Thus it would seem that PM codes are ideal for cosmo- 
logical N-Body simulations if we wish to study the large 
scale structure of the universe. 

The PM algorithm has been parallelised [s^] , typically 
this involves using parallel FFT 'si'l . The computational 
load is divided amongst processors by dividing the simu- 
lation box into slabs of equal size. 

Softening at the mesh scale ensures collisionless evolu- 
tion, but this also means that PM codes cannot resolve 
structure at length scales smaller than a mesh. Even in 
dense regions force softening is at the grid scale. This se- 
riously limits the effective dynamic range of simulations 
run with PM codes. A related problem is that the use of 
mesh makes force of a particle anisotropic at small scales. 

Improving force resolution in high density regions can 
improve the effectiveness of PM codes. Several tech- 
niques have been proposed for achieving this as poor 
resolution is the main shortcoming of PM codes for cos- 
mological N-Body simulations. 



5. Adaptive Mesh Refinement 

In Adaptive Mesh Refinement (AMR) the grid is re- 
fined in high density regions. A new mesh with smaller 
spacing is introduced and the low resolution force calcu- 
lated using the coarse global gri d is improved upon using 
the refined mesh 82, 83, 84, 85]. Several levels of refine- 
ment can be introduced; indeed are required in order to 
resolve substructure in dense haloes [SJ, |83 • Care is re- 
quired to ensure conservation of momentum and angular 
momentum in AMR codes. 



6. P^M: Particle-Particle + Particle-Mesh 

The basic idea here is to add a "correction" to the 
force computed using the PM method. This correction 
is computed by summing the contribution of close neigh- 
bours using the particle-particle method, hence the name 
PP-t-PM ~ P^M. It is assumed that this correction de- 
pends only on the distance, i.e., it is assumed to be 
isotropic, and is generally added out to a distance of 
about 1.5 — 2 times the distance between neighbouring 
mesh points^^ll^. P-^M was the first method to be used 
for high resolution cosmological N-Body simulations. 

The PP part of the calculation can also be done with 
the GRAPE chip |0 , though it requires some innovation 
as the chip is designed to return the force and not the 
short range correction. 

The P'^M has been parallelised, though load balanc- 
ing for such a code is not very simple as the overdense 
regions that require more CPU time are not distributed 
uniformly |88l l89ll90l| (see below). 

The P^M has some undesirable features. 



• The correction for the force is assumed to be 
isotropic, whereas the standard PM force has 
anisotropics at grid scale due to the anisotropic 
mesh structure. Thus the resulting force (long 
range -I- short range correction) must be anisotropic 
at the grid scale. 

• The short range correction in force is added only 
up to 1.5 — 2 grid lengths, whereas the PM method 
underestimates the force out to a much larger 
distance 1?^. 

• The refined inter-particle force is softened at scales 
much smaller than the average inter-particle sep- 
aration, this can lead to two body scattering and 
relaxation [ssIIb^ . Results of astrophysical interest 
like the correlation function may also get modified 
in the process 91] . 

• P'^M simulations slow down at late times when the 
distribution of particles becomes highly clustered. 
At this stage computation of the short range cor- 
rection of force dominates the total number of com- 
pute operations required. 

While the last item here relates to the efficiency of 
the code, other items in the list are far more serious as 
these raise doubts about the accuracy with which force is 
calculated. In order to retain the good features of P^M 
codes and address some of the issues listed above, several 
variations on the theme have been suggested. 

7. Tree + PM ^ TPM, GOTPM, TreePM . . . 

In this section we discuss a series of hybrid codes that 
combine the PM and the tree method in the same spirit 
as the PP and PM methods are combined in the P^M 
code. We will discuss these in order of the significance of 
departure from the P'^M method. 

• The Grid Of Trees PM (GOTPM) code [13 re- 
places the PP part of P^^M codes with a local tree in 
each region. This speeds up calculation of the short 
range force correction and the time taken for this 
calculation is less sensitive to the degree of clus- 
tering. This code takes care of the last undesirable 
feature listed for P'^M codes but it does not address 
any of the other issues. GOTPM is a parallel code 
where multiple levels of decomposition is used to 
achieve load balance [o^ . 

• The TPM code jUilliil also addresses the prob- 
lem of two body scattering in P'^M codes. In this 
code the short range correction to force is added 
only if the particle is in a high density region. Den- 
sity is computed at the position of each particle and 
a local tree is constructed in high density regions 
for computing the short range correction to the long 
range PM force. If the high density regions have a 
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large number of particles then the region is frag- 
mented into sub-regions. TPM is a parallel code, 
ideally suited for distributed parallel computing. 

• The Tree + PM (TreePM) code 96] makes sig- 
nificant changes to the P'^M approach. The force 
is partitioned between a short range and a long 
range component instead of using the force com- 
puted with a PM code as the long range force. In 
the TreePM code the long range force is truncated 
below a certain scale and the long range force 
becomes very small below this scale. In particular 
the long range force is small at grid scale if is 
chosen appropriately and therefore anisotropics in 
the long range force are also less important com- 
pared to the P^M, GOTPM and TPM codes. The 
short range force is computed using a global tree. 
The short range force has to be taken into account 
out to larger distances than in a typical P^M code. 
By tuning the model parameter Vs and a few other 
parameters in the code, it is possible to keep er- 
ror in force below 1% for most of the particles [o^ . 
The TreePM code is fairly simple to implement as 
the mathematical model of this method is well de- 
fined [93, . This code has been parallelised 
in a manner similar to that used for tree codes with 
some provision for the long range force calculation 
|lOCl| . The TreePM code solves all the other prob- 
lems of a P^M code but it does not address the 
issue of two body scattering and relaxation. 

• The Adaptive TreePM (ATreePM) code [lOl| ad- 
dresses the problem of two body scattering in 
TreePM code by using an adaptive softening length 
instead of a constant softening length. This can be 
thought of as an equivalent of AMR without using 
a refined grid. Local number density of particles is 
used to determine the softening length for particles. 
In order to ensure momentum conservation, force is 
symmetrised for particles that are closer than the 
softening length of either one of the two particles. 
Determination of local density and explicit sym- 
metrisation add a computational overhead. This is 
offset to some extent by the speedup at early times 
when the softening length is large for all the parti- 
cles. Using a hierarchy of time steps optimises this 
method further. 

The variety of techniques used for computing gravita- 
tional force discussed above is evidence of the work being 
done in developing better algorithms. It is also evident 
that full use has been made of parallel computing in order 
to achieve good performance. 

It is important to construct ways of comparing perfor- 
mance of different N-Body codes. We list some sugges- 
tions here. 

1. Dynamic range: The range of scales over which in- 
teraction force is computed reliably. There is rarely 
any problem at larger scales so the dynamic range 



is typically determined by the smallest separations 
over which the force errors are small. A useful unit 
for this scale is the average inter-particle separa- 
tion. We can fix this scale by requiring that at all 
larger scales error in force due to one particle be 
less than 1%. For good codes, this scale should 
coincide with the softening length. 

2. Trajectories of particles: The code should integrate 
the equation of motion in a reproducible manner 
and momentum should be conserved. The N-Body 
code should reproduce well known results about the 
correlation function and other statistical measures 
in the quasi-linear regime. Time step should be 
much smaller than the crossing time of particles in 
dense haloes, the ratio of time step to the smallest 
crossing time is a good estimate and this number 
should be smaller than 10~^. 

3. Efficiency: The N-Body code should be efficient 
and we should be able to run large simulations in 
as little time as possible. This requirement is likely 
to conflict with the first two requirements, and one 
should compare both the dynamic range and effi- 
ciency. It has been proposed that error in force 
should be plotted against time taken for comput- 
ing force for comparing codes we believe this 
to be the correct approach. 

4. Resource requirement: The N-Body code must 
store positions and velocities of all the particles, 
i.e., at least 6 numbers per particle. Most advanced 
codes discussed here store many more numbers per 
particle in order to speed up force calculation. This 
requirement can restrict our ability to run large 
simulations as the memory on computers is limited. 
In case of parallel codes, the relevant quantity is the 
maximum memory utilisation per particle for one 
processor. 

This ends the discussion of different algorithms for 
computing force in cosmological N-Body simulations. 



C. Setting up Initial Conditions 

N-Body simulations are generally started from fairly 
homogeneous initial conditions, i.e. the density contrast 
is much smaller than unity at all scales of interest in the 
simulation. In this regime we can use linear perturbation 
theory to compute all quantities of interest. In linear the- 
ory, the evolution of density contrast can be described as 
a combination of a growing and a decaying mode. At 
late times only the growing mode survives and hence we 
must choose the initial density and velocity field to put 
the system in the growing mode. Density contrast d is re- 
lated to 4> through eqn|21 and in linear theory the velocity 
field can also be expressed in terms of the gravitational 
potential in the growing mode. Therefore our problem 
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reduces to generating the gravitational potential and us- 
ing it to set up the density contrast and velocity field. 
The relation between the velocity and the potential in 
Zel'dovich approximation 18] is the same as that in the 
linear theory, therefore it is often said that the Zel'dovich 
approximation is use d to set up the initial conditions for 
N-Body simulations |102| |. 

Given an initial gravitational potential 0i„, there are 
two schemes for generating the initial density field. 

• The particles are distributed uniformly and their 
masses are chosen so that M = p(ti„)(l + Sin)dV, 
where Sin is evaluated at the position of the parti- 
cle. Here dV is the volume in the simulation box 
per particle and p(tin) is the average density at 
time tin- We can either start with zero velocities, 
in which case we have to increase the amplitude 
of (j)in by a factor 5/3 to account for the presence 
of decaying mode. Alternatively we can choose to 
put the system in the growing mode and assign the 
appropriate velocity to each particle. 

• Starting with a uniform distribution the particles 
are displaced using the velocity in linear theory 
for the growing mode. It is important to ensure 
that the maximum displacement is smaller than the 
average inter-particle separation in the simulation 
box. It can be shown that the resulting distribu- 
tion of particles will represent the required density 
field 123 if the initial distribution did not have any 
inhomogeneities. We can retain the initial velocity 
field used to displace particles. If the amplitude 
of displacements used is larger than the average 
inter-particle separation, it becomes necessary to 
recompute the potential from displaced positions 
and assign initial velocities with this potential (i^ . 
Such large displacements can lead to an incorrect 
realisation of the power spectrum. 

Schemes outlined above require an initial uniform dis- 
tribution of particles. This is important as any inhomo- 
geneities present in the initial distribution will combine 
with the density perturbations that are generated by dis- 
placing particles and will modify the initial conditions. 

• The commonly used solution is to place particles 
on a cubic grid, this is a uniform but not a random 
distribution. 

• An intuitive solution is to put particles at random 
inside the simulation box. This distribution has ^/n 
fluctuations which result in spurious clustering that 
tend to dominate over the fluctuations we wish to 
simulate. 

• Particles are placed in lattice cells but at a random 
displacement from the centre of the cell [t^ . This 
removes the regularity of grid without sacrificing 
uniformity. The amplitude of fluctuations can be 
controlled by reducing the amplitude of displace- 
ment about the centre of the cell. 



• The Glass initial conditions are obtained by evolv- 
ing an arbitrary distribution of particles in an N- 
Body simulation with a repulsive force. It can be 
shown that the amplitude of perturbations oscil- 
lates and decreases as a~^/^ [791. 



1. The Initial Gravitational Potential 

The initial density field is taken to be a Gaussian ran- 
dom field in most models. Linear evolution does not 
modify the statistics of density fields except for evolv- 
ing the amplitude of perturbations. As the potential and 
density contrast are related through a linear equation, it 
follows that the gravitational potential is also a Gaussian 
random field. A Gaussian random field is completely de- 
scribed in terms of its power spectrum ^103j . The Fourier 
components of a Gaussian random field (both the real 
and the imaginary part) are random numbers with a nor- 
mal distribution with variance proportional to the power 
spectrum of the random field. This property is used to 
generate the Gaussian random field in Fourier space and 
an inverse transform gives us the initial potential in real 
space. 

If some special features are required in the initial con- 
ditions, e.g., if we want a large planar perturbation, then 
we need to impose c onstraints on the G aussian random 
field to be generated [lOl [lOl [lOl IiOTI Gaussian ran- 
dom fields with a variable resolution |108| are needed for 
adaptive mesh refinement codes. 



III. LIMITATIONS OF N-BODY SIMULATIONS 

Here we discuss the scope and limitations of N-Body 
simulations. We will consider several issues concerning 
the domain of validity for N-Body simulations. 

N-Body simulations take initial density fluctuations 
over a finite range of scales into account. Do the fluctu- 
ations that are not taken into account make a difference 
to results of N-Body simulations? 

Several N-Body studies have shown that fluctuations 
at small sc ales do not affe ct structure that forms at large 
scales flOiilTTlllTTlITT^ in a signiflcant manner. Effects 
are of course there if larger scales are in linear regime 
but by the time larger scales reach 5^1, influence of 
smaller scales is ignorable. These studies used power 
spectrum, correlation function and visual appearance to 
reach this conclusion. Therefore in any N-Body simu- 
lations that are used, the scale where root mean square 
fluctuations are unity should be clearly resolved for re- 
sults to be reliable and independent of the small scale 
cutoff. On the other hand we expect some effect of small 
scal e fluctuations on how la rger density perturbations re- 
lax [nl [mini [mini, we can conclude that smaU 
scale fluctuations do not influence large scales in a sig- 
niflcant manner but it is an important issue and further 
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FIG. 1: Root mean square fluctuations in mass a are shown 
as a function of scale. The amphtude of g is for the ACDM 
model with = 0.7, Q.b = 0.05, h = 0.7 and n = 1.0. 
The linearly extrapolated amplitude is plotted here and no 
non-linear corrections have been included in generating this 
plot. 



studies are required to make the effect or the lack of it 
more quantitative. 

N-Body simulations assume that the density in the 
simulation box is same as the average density in the uni- 
verse. Therefore we must choose a simulation box such 
that the amphtude of fluctuations in the universe (or 
the model being simulated) at that scale is ignorable. 
Studies have shown that violating this requirement leads 
to an underestimate of correlation function though the 
mass f unct i on o f small mass haloes does not change by 
much Ill9l ll2Cl . Effects at large scales can be signifi- 
cant jl2ll 1123 1!^ In other words, the formation of small 
haloes is not disturbed but their distribution is affected 
by non-inclusion of long wave modes. One way of quan- 
tifying fluctuations at large scales, and hence their effect 
on structure formation at small scales is the amplitude of 
fluctuations at the scale of the simulation box. Figure [21 
shows lines of constant root mean square fluctuations in 
mass a. These lines are plotted as a function of scale 
and redshift for the ACDM model (f^A = 0.7, = 0.05, 
h = 0.7 and n = 1.0). Curves are for a = 0.1, 0.05, 0.025 
and 0.01 in ascending order. We can flnd out the lowest 
redshift up to which a given simulation box may be used 
once we flx a that is considered acceptable at the box 
size. 

However, a threshold in a does not carry any infor- 
mation about the shape of the power spectrum and that 
is extremely relevant here. We find that the best way 
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FIG. 2: Lines of constant root mean square fluctuations in 
mass a are plotted as a function of scale and redshift for the 
ACDM model (see text for details). Curves are for a = 0.1, 
0.05, 0.025 and 0.01 in ascending order. We can find out 
the lowest redshift up to which a given simulation box may 
be used once we fix a that is considered acceptable at the 
box size. Such a threshold does not carry any information 
about the shape of the power spectrum. Comparing with an 
approach based on convergence of collapsed mass, we fi nd that 
(7 = 0.025 is a reasonable choice at low redshifts |123| . 

of quantifying the effect of long wave modes is to check 
whether including them in the simulation will change the 
number of massive haloes or not I2.''> . This can be es - 
timated using the Press-Schechter mass function |l24j . 
If there is a measurable effect on the number of haloes 
and collapsed mass, the long wave modes are relevant 
and must be considered. The large scale structure in N- 
Body simulations does not c onverge until all such modes 
are taken into account |12,'^ . Our results for the ACDM 
model are: 

• A box size of 150h~^Mpc is needed for simulations 
that are evolved to the present epoch {z = Q). 

• Simulations that are evolved up to z '--^ 3 should be 
50h^^Mpc across for generic applications. 

• Simulations for studying early reionisation (z ~ 10) 
should be at least 20h~^Mpc across if these are to 
have a representative density field. 

Comparing these with the curves in figure[21 we find that 
a — 0.025 is a reasonable choice except at very high 
redshifts. a = 0.01 is a better choice at higher redshifts, 
and this variation in the threshold value is indicat ive o f 
the dependence on the shape of the power spectrum jl23j | . 
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The effect of a finite simulation box on modes com- 
parable with the box size can be estimated analytically 
|l25j and this can also be used to check whether the cho- 
sen physical size of the simulation box is sufficiently large 
or not. 

Methods have been de veloped to take the missing long 
wave modes into account |l26lll27j . these are wave modes 
larger than the simulation box. These methods evolve a 
realisation of long wave modes independently and their 
effect is combined with the evolution of small scales in the 
N-Body simulation. This allows one to work with smaller 
simulation boxes, but the methods have limitations and 
cannot be used for simulating very small regions. Indeed, 
the criterion based on convergence of collapsed mass can 
be used to fix the smallest box size for which these meth- 
ods can be employed (12,'^] . 

Other limitations relate to two body relaxation, poor 
accuracy in force calculation and discreteness. It has 
been asserted that in coUisionless simulations forces are 
not required to be very accurate as the di screteness noise 
dominates over errors in force 12^ Il29j . This is more 
relevant for non-cosmological simulations, where many of 
these tests were carried out. More recent studies of force 
softening have pointed out that the form of softening 
can be very important [s^ in controlling the discreteness 
noise. 

Tests that are specific for cosmological simulations |9ll | 
show that unless mass and force resolution are compa- 
rable, two body relaxation can leave detectable signa- 
tures. Thus it is important to avoid discreteness effects 
and maintain parity between mass and force resolution 
in N-Body codes. 



IV. NON-LINEAR GRAVITATIONAL 
CLUSTERING 

N-Body simulations have been used extensively to un- 
derstand aspects of gravitational clustering. Of the four 
fundamental interactions, gravity is the weakest and it is 
impossible to carry out laboratory experiments in grav- 
itational dynamics of a many body system as other in- 
teractions overwhelm gravity at these scales. N-Body 
simulations can be used as a test-bed for doing numer- 
ical experiments in order to understand various aspects 
of non-linear gravitational clustering. Indeed, we cannot 
claim to have understood the process of structure for- 
mation till we develop sufficient insight so as to make 
N-Body simulations redundant except for the purpose of 
computing detailed predictions for models of structures 
formation. 

Considerable progress has been made using quasi- 
linear approximations, non- linear scaling relations and 
N-Body simulations, though many important questions 
are yet to be settled. Basic premise here is that if we are 
interested in only a limited amount of information about 
the final state then it should be possible to simplify the 
calculation by making an ansatz that captures the essen- 



tial physical process at work. One can also invert this 
and find out which physical process dominates in a given 
situation. We review some of the issues that have been 
studied in detail. 

The distribution of matter is very close to homogeneous 
at early times, how is it confined to central region of 
potential wells at late times? 

• Infall is described very well by first order La- 
grangian perturbation theory. This, extrapolated 
to mildly non-linear density contrasts is called the 
Zel'dovich approximation [13 . This approximation 
suggests that the first structures to form as a re- 
sults of collapse are typically surfaces of high den- 
sity, the so called pancakes. Comparisons with N- 
Body simulations show this t o be a very good ap- 
proximation up to this stage |l30l |. However, these 
surfaces of high density thicken and disappear if 
we extrapolate the approximation scheme to late 
times, whereas these pancakes do not thicken in N- 
Body simulations. Clearly, the Zel'dovich approxi- 
mation and even higher ord er Lagrangian perturba- 
tion theory |l3ll Il32l Il33j | lacks the key ingredient 
that helps confine matter to potential wells. 

• An artificial viscosity term can be added to the 
equation of motion, eqn.l^, and it can be trans- 
formed into Burger's equation p^ . This prevents 
interpenetration of infalling streams and pancakes 
remain thin, this is called the Adhesion approxi- 
mation |l9l |20| . The approximation puts matter in 
the right place but fails to find a physical reason 
for the viscosity term. 

• The equation of motion, eqn.JQ), contains a velocity 
dependent term where the expansion of the uni- 
verse acts as a viscous drag. Further, it can be 
shown that the gravitational potential varies very 
slowly. If we assume that the potential changes 
at the rate expected in the linear theory while re- 
taining its initial spatial dependence then the drag 
force is sufficient to confine matter to the central 
regions of potential wells J2, 23] . This approxi- 
mation {Linear Evolution of Potential (LEP), aka 
Frozen Potential Approximation (FPA)) compares 
well with N-Body simulations. We can conclude 
that the drag term in the equation of motion plays 
an important role in confining matter to potential 
wells and the interaction of infalling particles plays 
a less significant role. 

How strongly do density fluctuations at small scales in- 
fluence density perturbations at large scales? This is 
relevant as N-Body simulations ignore perturbations at 
scales smaller than those probed in the simulation. 

• If there are no perturbations at large scales then 
small scale fluctuations generate a P{k) oc k* spec- 
trum at small k, this effect can be derived in the 
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second order perturbation theory as well as from 
the full set of equations . Amplitude of this limit- 
ing spectrum increases rapidly at early times |lllj . 
This effect of course cannot be seen if there are fluc- 
tuations at large scales with P{k) oc fc" with n < 4 
at small k. 

• If there are small density perturbations present at 
large scales then there is a discernable visual ef- 
fect of small scales o n pe rturbations at large scales 
[lOgl UTol [ml ITTl In^ for generic initial condi- 
tions. Small scales can influence the power spec- 
trum and higher moments of particle distribution 
at these scales, the detailed effect depen ds on the 
form of fluctuations at small scales |l34j . 

• If there are non-linear density perturbations 
present at large scales then there is no discern- 
able visual effect of small scales o n pe rturbations 
at large scales 'lol [ll3, lUll Ell El for generic 
initial conditions. There is no significant effect in 
the powe r spectrum or the two point correlation 
function |l09l llllj in such a case. If the small 
scales and the scale of non-linearity are well sep- 
arated then there is no effect in higher moments 
either p^ . 

• Clearly, non- linear gravitational clustering works to 
remove the effect of perturbations at large scales. 

• For special initial conditions, it has been found that 
collisions between clumps formed due to small scale 
perturbations can lead to a faster relaxation of per- 
turbations at large scales. E.g., if larger scales are 
modelled as a single plane wave then the result- 
ing pancake is thinner in p resence of significant 
small scale fluctuations |118| . It remains to be seen 
whether this is an important effect for generic ini- 
tial conditions. 

Does gravitational clustering erase memory of initial con- 
ditions? Is there a one to one correspondence between 
some characterisation of initial perturbations and the fi- 
nal state? Note that this is different from removing ef- 
fects of perturbations at scales much smaller than the 
scale of non-linearity. Clearly, if the answer to the first 
question is yes then we cannot recover any information 
about the initial density perturbations from determina- 
tion of density perturbations in the non-linear regime. 

• N-Body simulations show that gravitational clus- 
terin g do es not erase memory of initial conditions 
[llllillllsllIlllIlP,'!^- The final power spec- 
trum is a function of the initial power spectrum and 
this relationship can be written as a one step map- 
ping. 

• The functional form of thi s ma pping depends on 
the initial power spectrum |140| . 



• An analytical understanding of some aspects of this 
mapping can be deve l oped using simple models and 
approximations |l4llll42| |. 

• It is found that density profiles of massive h aloes 
have a fo rm independent of initial conditions |l43l 
Il44l Il45| . It is important to note that there is 
considerable scatter in density profiles obtained 
from N-Body simulations and it is difficult to state 
whether a given functional form is always the best 
fit or not. T here are a large number of recent stud- 
ies (e.g. see fT4^IT47| 'l but most of these focus on 
the ACDM model and do not explore other initial 
conditions. 

Is it possible to predict the masses and distribution of 
haloes that form as a result of gravitational clustering? 

• The initial density field is taken to be a Gaussian 
random field, and for hierarchical models the 
simple assumption that each peak undergoes col- 
lapse independent of the surrounding density dis- 
tribution can be used to estimate the mass function 
fl24, 148] and several related quantities. 

• N-Body simulations show that this simple set of 
approximations is incorrect, however the resulting 
mass function is fairly accurate over a wide range 
of masses. 

• Merger rates can be computed using the extended 
Press- Schechter for malism Il48l. t hese are useful for 
many applications |l49l 115(1 Il5l| . It is also possi- 
ble to estimate clustering properties of haloes using 
this formalism |l52lll53| . 

• Modifying some of these assumptions can lead to 
improved predictions |l54ll55l.ll56l |. 

We have highlighted some of the key issues in non- 
linear gravitational clustering in this section and re- 
viewed how N-Body simulations and various approxima- 
tions have been used to develop insight. N-Body simula- 
tions have been and are being used for a variety of other 
applications, the most notable being computing detailed 
predictions for models of structure formation. 

V. GASTROPHYSICAL EFFECTS 

N-Body simulations that take only gravitational in- 
teractions into account can be used to obtain the large 
scale distribution of galaxies and clusters. Further details 
can only be obtained using simulations with gas physics. 
There are two very different approaches to including gas 
physics in N-Body simulations. 

The classical approach is to use a gri d to solve fluid 
equations \l57, 158, 151 [US [HI UHl fl64, 165i 
Il67l Il68l| . Fluid interactions are short range and infor- 
mation from nearby grid points is sufficient to evolve fluid 
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properties at any grid point. Of course, fluids must re- 
spond to gravitational field of the matter distribution. 
This type of a code is well adapted for capturing shocks 
and discontinuities. Resolution of the code can be en- 
hanced in high density regions by using mesh refinement. 
In cosmological simulations it is important to improve 
the mass resolution of dark matter particles along with 
an enhancement in resolution for hydrodynamics. Grid 
codes can be easily generalised to include other effects 
such as magnetohydrodynamics |l69l llTOj . Such codes 
have been used effectiv ely fo r cosmological applications 
[T57lll6lll6llTMIlMll67llT68j . 

The Smooth Particle Hydrodynamics (SPH) pTllIlT^ 
where particles are assigned fluid properties is a very dif- 
ferent approach to the same problem. The fluid prop- 
erties like pressure, density, temperature, etc. at any 
point can be found by averaging over particles in the re- 
gion using a weight function. Use of interpolation for 
determining fluid properties makes it impossible to re- 
solve shocks and discontinuitie s in SPH codes. There 
are several known limitations |l73l Il74l Il75| that one 
should be aware of. SPH codes are relatively easy to im- 
plement and several implementations are in use. Several 
implementatio ns for cosmological simulations are in use 
1 1781. Va riations of SPH with focus on entropy 
equation 17y, or on simulation of multi-phase media 
[18I,, ,182J have been developed. 

Both typ es of codes have been compared and give sim- 
ilar results |l83l Il84 l| for astrophysical applications. 

Effects other than hydrodynamics can also be incor- 
porated in a similar fashion in both types of codes. For 
example elementary chemical reactions like formation of 
Hydrogen molecules, cooling, heating, etc. are impor- 
tant local effects. Key effects like star formation and 
feedback from stellar and other sources are difficult to 
include as vastly different scales are of relevance. As a re- 
sult, much of the treatment of these effects has remained 
phenomenological. 

Radiation transpor t is a non-local effec t and is di f ficul t 
to take into account jTH [HI [HI [HI [TH [l9i [ip. 
The radiation field is either assumed to be isotropic and 
homogeneous, or it is assumed to originate at a few point 
sources. The time scales over which the radiation field 
evolves are much shorter than most other time scales in 
the problem. Thus we can assume that the density and 
velocity field do not change while studying evolution of 
the radiation field. 

State of the art simulations can be used for studying a 
variety of problems. We summarise three systems here: 



• The inter-galactic medium [T9i [M [T9l[T95l[iM 
Il97j is believed to contain mostly pristine mate- 
rial with elements other than Hydrogen and Helium 
present in only trace amounts. The radiation field 
can be assumed to be isotropic and homogeneous at 
z < 3, though at higher redshifts it is patchy. Den- 
sities in the inter-galactic medium do not become 
very large. 

• The intra-cluster medium IM ^84, 198', [TH [203 
I201I I20I 12ml I20I l2nl ma [207.1 is believed to 
be close to hydro-static equilibrium. Gas in the 
intra-cluster medium is so hot that many radiative 
processes can be ignored safely. There are several 
sources of energy in clusters of galaxies, most of 
these are distributed uniformly with some contri- 
bution from a central AGN. Magnetic fields may 
be important enough to make the physics of the 
problem more complicated. 

• Formation of first stars [20I [20I [2l3, [ml 

[21I [2TI [2TI requires simulating cooling 

of Hydrogen (mainly) through various processes 
till molecules form and lead to temperatures low 
enough for formation of stars. 

Focus in future will be on incorporating further effects 
so that more complex problems like galaxy formation 
can be studied in full detail. N-Body simulations are 
useful not only as tools for evolving complex systems, 
these can also be used to understand which effects play a 
more important role in different phases of this evolution. 
Next decade holds the promise of understanding physics 
of galaxy formation with help of N-Body simulations. 
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